Equilibrium and stability of supermassive stars in binary systems 
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We investigate the equilibrium and stability of supermassive stars of mass M > 1O 5 M0 in binary 
systems. We find that corotating binaries are secularly unstable for close, circular orbits with 
r < AR{M yiO 6 Mq) 1 / 6 where r is the orbital separation and R the stellar radius. We also show that 
corotation cannot be achieved for distant orbits with r > 12R(M /10 6 M©) -11 / 24 , since the timescale 
for viscous angular momentum transfer associated with tidal torques is longer than the evolution 
timescale due to emission of thermal radiation. These facts suggest that the allowed mass range and 
orbital separation for corotating supermassive binary stars is severely restricted. In particular, for 
supermassive binary stars of large mass M > 6 x 10 6 Mq, corotation cannot be achieved, as viscosity 
is not adequate to mediate the transfer between orbital and spin angular momentum. One possible 
outcome for binary supermassive stars is the onset of quasi-radial, relativistic instability which drives 
each star to collapse prior to merger: We discuss alternative outcomes of collapse and possible spin 
states of the resulting black holes. We estimate the frequency and amplitude of gravitational waves 
emitted during several inspiral and collapse scenarios. 



I. INTRODUCTION 



Recent astronomical observations provide increasingly strong evidence that supermassive black holes of mass ~ 
10 6 — 10 10 M Q , exist and that they are the central engines of active galactic nuclei and quasars (l). In addition to 
their importance for a fundamental understanding of active galactic nuclei and quasars, gravitational waves from the 
formation of the supermassive black holes and from inspiraling supermassive binary black holes are likely sources for 
proposed laser interferometric gravitational wave detectors in space, such as the Laser Interferometer Space Antenna 
(LISA) J2|,[|. However, the scenario by which supermassive black holes form is still uncertain. Viable stellar dynamical 
and hydrodynamical routes leading to the formation of supermassive black holes have been proposed (see, e.g., [0). 
In one hydrodynamical scenario, a supermassive gas cloud is build up from the multiple collisions of stars or small 
gas clouds in stellar clusters to form a supermassive star. Supermassive stars subsequently collapse to black holes 
after quasi-stationary contraction as a result of cooling |^]|3). Such a scenario for forming a supermassive star and 
black hole has been investigated in detail in the 1960s and 70s using simple analytical models or assuming spherical 
symmetry (see, e.g., for a review and references). 

Supermassive stars are likely to be rapidly rotating (see, e.g., [§[). For studies of rapidly rotating stars in equilibrium, 
it is necessary to perform numerical computations using modern numerical techniques j^] which were not feasible 
until quite recently. Baumgarte and Shapiro recently performed a detailed numerical analysis of a rapidly rotating 
supermassive star J7J. Because of the fact that the viscous or magnetic braking timescale for angular momentum 
transfer is likely to be shorter than the evolution timescale of supermassive stars in the typical mass range, the star 
can be assumed to be in rigid rotation at the mass-shedding limit (which implies that the spin angular velocity of each 
star is equal to the Kepler velocity) JTo) ] . They find that at the onset of quasi-radial collapse, the non-dimensional spin 
parameter cS/GM 2 , where S,M,c and G are spin, mass, light velocity, and gravitational constant, is nearly equal to 
unity, so that rotating supermassive stars might not collapse directly to a black hole. They speculate that the collapse 
may be inhibited, forming a disk and/or a bar, or possibly fragmenting into several blobs. If such non-axisymmetric 
structures form during the collapse, quasi-periodic gravitational waves of ~ 10~ 4 — 10~ 2 Hz may be emitted and 
provide a strong source for laser interferometeric gravitational wave detectors in space Q . 

In this paper, we consider supermassive stars in binary systems. Such binaries may arise naturally from fragmenta- 
tion of a large gas cloud with rapid rotation, or during multiple mergers of small gas clouds and/or stars [[4|]S|]. Even 
in binary systems, each supermassive star is likely to reach quasi-radial instability eventually after quasi-stationary 
contraction. In binaries, most of the angular momentum of the system can be distributed in orbital motion so that 
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the spin of each star may be much smaller than that for isolated supermassive stars. In particular, small spin can be 
achieved for the case when viscous transport between orbital and spin angular momentum is effective. If so, each star 
is in slow rotation with respect to the inertial frame, and collapse to a black hole may proceed in an almost spherical 
manner, with tiny emission of gravitational waves. Thus, the scenario for forming black holes in binary systems could 
be considerably different from that for isolated supermassive stars. Also, a binary containing supermassive stars can 
be a progenitor of an inspiraling, supermassive black hole bin ary in a close orbit. Such a binary can emit gravitational 
waves in the frequency range appropriate for LISA [see Eq. ( |3.l| )]. To detect gravitational waves from binaries using 
matched filtering technique ||ll| , |l2f , it is useful to anticipate possible mass and spin states of the black holes. With 
this motivation, we investigate the equilibrium and stability of supermassive binary stars in this paper. 

The paper is organized as follows. In Sec. II, we estimate various timescales which are relevant for evolution of 
supermassive binary stars. In Sec. Ill, we construct analytical equilibrium models of corotating, supermassive binary 
stars and investigate their stability. In Sec. IV, we present numerical solutions for corotating, supermassive binary 
stars to demonstrate that the analytical results shown in Sec. Ill are qualitatively correct. In Sec. V, we consider the 
evolution of corotating, supermassive binary stars as a result of thermal radiation, and clarify their final fate. In Sec. 
VI, we briefly discuss gravitational waves from supermassive stars in the present scenarios. Section VII is devoted to 
a brief summary. Hereafter, we adopt the geometrical unit in which G = c = 1. 



II. TIMESCALES 



In this paper, we investigate equilibrium and stability of supermassive stars in close binary systems. Before proceed- 
ing, we estimate the relevant timescales characterizing such systems to clarify possible states of supermassive binaries. 
Here, we define supermassive stars to be self-gravitating objects of mass larger than ~ 10 5 M Q (see discussion below). 
Also, we assume that the binary is an isolated system; namely, the accretion of mass and heating from outside are 
assumed to be negligible. 

The equation of state for supermassive stars is dominated by thermal radiation pressure with a small correction 
due to baryonic gas pressure. Such stars are convective with constant entropy per baryon |^|,^|J^]. These conditions 
imply that the structure of the supermassive stars is very close to that of an n = 3 polytrope so that 

P = Kp^ 3 + P g , Pg^Kp 4 / 3 , (2.1) 

where P, p, P g , and K are the total pressure, baryon mass density, gas pressure, and polytropic constant, respectively; 
K is approximately related to the entropy of radiation, s r , according to ||U 
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Here, a is the radiation density constant, and mg the mass of hydrogen; we assume that supermassive stars are 
composed of pure hydrogen for simplicity. For n = 3 spherical polytropes in Newtonian theory, the mass M is 
uniquely determined by K according to 



m = (ttJ ■ <"> 

where k\ and &2 are constants of order unity associated with the Lane-Emden function (see Sec. III). Using this 
relation, s r can be written approximately as H 

/ M n 1/2 



where k denotes the Boltzmann constant. 

In Newtonian theory, pure n = 3 spherical polytropes are marginally stable to radial collapse, and with a tiny general 
relativistic (GR) correction, they are unconditionally unstable. To stabilize such a star, some additional perturbations 
are necessary. One perturbation is spin angular momentum. As shown in [Q, supermassive stars are stabilized if they 
have a large spin angular velocity. However, in this paper, we mainly analyze corotating, supermassive binary stars, 
whereby the spin angular momentum of each star is not large enough to stabilize the star against radial collapse. To 
stabilize such slowly rotating supermassive stars, the role of baryonic gas pressure, which constitutes a small fraction 
of the total pressure, is crucial. 
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The evolution of isolated supermassive stars that do not have large angular momentum and are stabilized by the 
gas pressure has been investigated in detail many years ago |]||. According to a typical scenario, a supermassive 
star shrinks quasi-statically as a result of emission of the thermal radiation from its surface until it reaches a critical 
point where a GR instability sets in against radial collapse. The radial instability occurs for a critical compactness 
IH [see also Eq. j ^2% ] 

-> 0.445 (T--y (2.5) 

where R is the circumferential radius, and T is the effective adiabatic constant. Here, T may be approximately written 
as 

where (5 is the ratio of the gas pressure to the radiation pressure, (3 = P g / (Kp 4 / 3 ) <C 1, which can be expressed as || 

Thus, the constraint of M/R for a stable supermassive star of mass M can be rewritten as 

M ( M \ ~ 1/2 

« S0JS W ■ < 2 - 8 ' 

In this evolutionary scenario, there are two relevant timescales; the emission timescale of thermal radiation and 
the dynamical timescale. Supermassive stars emit the thermal radiation from the stellar surface with the Eddington 
luminosity 

r 4itM I* c\\ 

^Edd = , (2.9) 

where kt is the opacity with respect to the Thomson scattering. The total energy of a supermassive star is approxi- 
mately §§ 

M 2 

(2.10) 

Thus, the emission timescale of thermal radiation is 
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E\ , „ inl0 (M/R\f 



-^Edd 

where, in evaluating kt, we assume pure ionized hydrogen. The dynamical timescale is approximately 



rd yn .^/^1.4xl0 4 sec(^) ^(^l, (2-12) 



V 1O 6 M 

where p c denotes the maximum (central) density of supermassive stars which, for n — 3 spherical polytropes, can be 
expressed as p c ~ 12.94M/R 3 (see next section and Table I). 

For a star to evolve quasi-statically, the inequality, t cvo i > Tdyn, must be satisfied, so that 
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where we use Eq. (2/7) for evaluating /3. From the requirement of stability, M/R also has to satisfy Eq. (|2.q). Thus, 
for existence of a stable supermassive star in quasi-static evolution, the mass has to be less than ~ 2 x 10 8 M Q . If 
the mass is larger than this value, supermassive stars collapse dynamically without evolving through a quasi-static 
phase. Furthermore, supermassive stars of mass M < lO 5 Af0 reach sufficiently high temperature for nuclear burning 
to become important before reaching the onset of the instability. If nuclear reactions set in, the temperature increases 



3 



so rapidly that the star undergoes a violent explosion jL4|. Thus, we restrict our attention to stars with M > 10 5 M© 
which are stable until they reach the onset of radial instability. 

For binary systems, several additional timescales have to be taken into account. One is the radiation-reaction 
timescale of gravitational waves, which is written using the quadrupole formula as [ [l5[ 

5 r 4 u fr/M\ A f M \ , 

^5i2M^ 2x10 sec Tor (10^ h (2 - 14) 



where we assume that the mass of each star, M, is identical (i.e., the total mass is 2M), and where r denotes the orbital 
separation between two centers of mass. Since r/M > R/M > 10 3 for stable supermassive stars with M > 10 5 M Q , 
we immediately find tqw 3> T dyn and tcw 3> T cvoh so that gravitational radiation reaction is not relevant in this 
problem. 

According to standard scenarios, supermassive stars were formed near galaxy centers in the early stage of galaxy 
formation Q. Associated with galaxy formation in its early stage, gas may settle toward the center of the galaxies. If 
the gas accretes onto supermassive binary stars, it can be a source of friction which can decelerate the orbital motion 
of the stars. The dragging timescale can be written crudely as 

M 4% f p c \ 1/2 fp m \ f r^ l/2 



Tdras ^~iUJ \TJ\r) Tdyn ' (2 - 15) 

where p and p m denote the typical density of the accreting gas, and the mean density of supermassive stars. If 
Tdrag < Tcvoi the orbital decay due to friction can become effective in inducing merging. However, to achieve this 
condition, po has to be extraordinary large, 



implying that po has to be larger than ~ 10 13 M Q /pc 3 . Also, the total amount of mass of the accreted gas in time 
r cvo i has to be as large as ~ M so that the binary in such situation cannot be regarded as an isolated system contrary 
to our assumption in Sec. II. Thus, we do not consider the effect of the friction. 

The viscous timescale plays an important role in determining the velocity field in binary systems. Here, we have to 
take into account two timescales. One is associated with the tidal torque from the companion star (hereafter, we refer 
to this viscous timescale as r v ; s ). If r V j S is shorter than r ovo i, transfer between orbital and spin angular momentum 
can be efficient and corotation is achieved. On the other hand, for t v ; s > r ovo i, angular momentum transfer is not 
important, implying that the spin of each star is determined solely by the viscous angular momentum transfer inside 
each star (hereafter, we refer to this timescale as T v j s o). In the case when t v ; s > r cvo i > t v ; s o, each star will be in rigid 
rotation in its comoving fr ame. For T V j S > r V i S o > t ovo i, vorticity of stars is conserved, but such a high value of r v i s o 
is unlikely [see below Eq. ( |2.24| )]. 

The timescale for the viscous transfer of the angular momentum induced by the tidal torque is approximately |lq] 
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(2.17) 



where v is the viscous parameter. As often found in astrophysical contexts, microscopic, molecular viscosity yields 
timescales that are larger than the evolutionary timescale by many orders of magnitude. However, the effect of 
turbulent viscosity can be estimated by assuming that the velocity of turbulent motion vt is an appreciable fraction 
of the velocity of sound JT^] , 

«t = CWsoimd, ( 2 - 18 ) 

where we take the dimensionless viscosity parameter a to be in the range 0.01 < a ^ 1. Assuming the characteristic 
length scale of turbulence to be of O(R), we can estimate v as 

v ~ aRv sound . (2-19) 

Then, 
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To achieve corotation, the condition t v - is < t cvo \ has to apply, and hence, 

5 / 12 / i\,r \ -V4 



Since M/R has to satisfy Eq. (J2.8|), we obtain 



4 <12cW- 
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(2.21) 



(2.22) 



R ~ V 106M ©. 
The existence of binary systems requires r to be larger than ~ 2R, implying 

M < 5 x lO 7 a 4/11 M . (2.23) 



Therefore, supermassive stars in corotating binaries can exist only for close orbits [cf. Eq. (2.21)], for compact stars 
[cf. Eq. ( p. 22 )], and for a fairly small mass [cf. Eq.(2.23)]. In particular, supermassive stars in binary systems with 



masses larger than ~ 1O 8 M0 cannot be in corotation. In the next section, we will show that the stability may further 
restrict the allowed region for these parameters. 

Even if tidal interactions are not important, the viscous transfer of angular momentum inside each star can be 
important j?J, since the viscous timescale is 

TvisO ~ a _1 Tdy n . (2.24) 
Thus, as long as M/R > 1.7 x lO~ 6 or 2/5 (M/lO 6 Af ) 3/5 , the relation T vis0 < t cvo i always holds. 

III. ANALYTICAL MODEL FOR COROTATING BINARY 

We investigate the stability of supermassive binary stars supported by a polytropic equation of state 

P = Kp T . (3.1) 

In this section, V is chosen to be close to 4/3 according to 

4 
r=-+e, e =£« L ( 3 - 2 ) 

With this choice of equation of state, we can investigate the quantitative properties of supermassive stars of mass 
M 3> Mq analytically by means of a variational principle. As mentioned in Sec. II, such stars have an effective 
adiabatic constant given approximately by 
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(3.3) 



If we assume that the mass of supermassive star s is in the range between 1O 5 M and 2 x 1O 8 M , e should be in the 
range between ~ 1CT 4 and - 5 x 1CT 3 [Eq. (O)]. 



We here consider corotating supermassive stars of equal mass for simplicity. The alternative extreme case of binary 
component with large mass ratios is also investigated in Appendix A using a Roche binary model. There we show 
that the results are qualitatively the same. To determine equilibrium and stability of supermassive stars in corotating 
binary systems, we write the total energy for each star as the sum of the internal energy U, the potential energy W, 
the orbital kinetic energy T, the spin kinetic energy T s and the binding energy between two stars Wb as H|l8f 

U = k 1 KMx n l' i , (3.4) 
W = ~k 2 M 5/3 x - k 4 M 7/3 x 2 , (3.5) 

T=fr 2 tf, (3.6) 

T s = \m 2 , (3.7) 
M 2 

W b = - — . (3.8) 
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Here, M is the mass of each star, x = p c the density parameter (p c is the maximum density again), r the orbital 
separation between the two centers of mass, n[= l/(T— 1)] the polytropic index, O the orbital angular velocity, and I 
the moment of inertia. We do not distinguish the rest mass (baryon mass) and gravitational mass since their difference 
is quite small and unimportant for supermassive stars. The nondimcnsional structure constants k\, k 2l and k^ are 
constants dependent on n calculated from the Lane-Emden function 11. In Table I, we list these constants for several 
n near 3. 

In W, we include the first post-Newtonian correction taking into account GR effects on the stability against quasi- 
radial collapse. We neglect the first post-Newtonian corrections in Wb since their magnitude is smaller than that of 
the Newtonian terms by several orders of magnitude as M/r < M/2R <C 1. We also neglect the tidal energy in Wb 
and W |19[ ). As we estimate the magnitude of this effect in Appendix B, it is not substantial for corotating binaries. 
As long as e > 7 x 1CP 4 , it is not important even in non-spinning binaries. 

The quantity I may be written in the form 

I = kM 5 / 3 P - 2 / 3 , (3.9) 

where k is a constant calculated from the Lane-Emden function (cf. Table I). 
The total energy per one star, E, then becomes 

E = k x KMx*' n - k 2 M 5 ' 3 x - k 4 M 7 / 3 x 2 + + n 2 - (3.10) 

The equilibrium angular velocity f2 is derived from the relation dE/dr\M,j,K — as 



n = J 2 ™ (3.ii) 



where J is the angular momentum of one star 



J = +m= (ryfT 1 / 3 + KX- 2 n)M 5 /\ (3.12) 

and r\ = 2~ 4 / 3 . We note that the second term denotes the spin of each star S(— 7X2). 

Substituting r = (2M/Q 2 ) 1 / 3 in E and taking the first derivative of E(= E[x,Q(x)]) with respect to x, fixing M, 
J and K, yields a second condition for equilibrium, 



ox 



= —^KMx 3e - k 2 M b ' 3 - 2fc 4 M 7/3 z + '\ (3.13) 

m,j,k n x 6 



where we use e = 1/n — 1/3 and the relation 



dx 



M.J X° \ 3 X 

The second derivative of E(= E[x, £l(x)]) with respect to x becomes 
d 2 E 9efa 



= ^fif_^-4/3 + «A 1 <Q (314) 



dx 2 



EMx 3 ^ 1 - 2k 
m,j,k n 



" x 1 \r]x 2 -3/d7 4 / 3 ) 



, x 7/ , 3k« 2 M 5 / 3 / rix 2 + k^ 3 \ 
= 3ek 2 M^ - 2(1 - 3e)fc4 MV3 __ _ 3^/3 + <) > W 

where to derive the second line, we use Eq. (E3T3). For stable equilibrium, this second derivative has to be positive, 
and the zero point marks the onset of instability ]|]6j . 

It is convenient to analyze equilibrium and stability fixing Vl 2 j p c = fl 2 /x 3 = lu 2 , which is equivalent to fixing the 
ratio of the orbital separation to the stellar radius of stars: 



2 Sir ( R\ ( p c 



= | S J ~,M»(22) , (3.16, 
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where p m denotes the mean stellar density, i.e., 3M/(47r_R 3 ), and p c j ' p m ~ 54 for e <C 1 (see Table I). Since r > 2R 
for binary systems, uj 2 should be less than ~ 0.02. (In reality, ui 2 is less than ~ 0.01 since the equator of each star is 
significantly enlarged by the centrifugal force associated with the spin for close binaries.) 
The equilibrium condition dE / 8x\m,j.k — yields the relation between M and x as 



M = 



-Qm - kuj 2 ) + y/(k 2 - KUJ 2 ) 2 + 24k 1 k 4 Kx 1 +^/r 



Ak 4 x 
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(3.17) 



The condition for the onset of instability d 2 E / '8x 2 \m,j,k — with Eq. ( 3.17 ) gives the solution of x at the critical 
points as 



3(efc 2 - Kcu 2 F a ){k 2 - klu 2 (1 + 3F - 3e)}l " /3 



2(l-3e) 2 (l + 3e)kik 4 K 



— Xc > 



where 



r\ — 3koj 4 / 3 



(3.18) 



(3.19) 



Since o ur pu rpose here is to determine the stability, we restrict our attention to the case where K is constant. 
Equation (3.18) implies that on every equilibrium curve, M(x), of a fixed value of lo, there exists one turning point 
x c for a sequence along which J is constant as 



J = (r)uj- 1/3 + klo)x- l/2 M(x c ) 



1/2; 
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(3.20) 



However, the implications of the existence of turning points for small u> and large u> are different. Consider the relation 
between M = MK~ n l 2 and (p c )( 3 - n )/ 2 " (p c = p c K n ) for e = 0.003 in Fig. 1 as an example. We remark that the 
qualitative behavior is independent of e as long as < e < 0.01. In the presence of a general relativistic effect or in 
the case e > 0, MK~ n / 2 is not constant even in the infinite separation (i.e., uj 2 — 0), but it is still close to ~ 4.5 
for a wide range of the central density and orbital separation as long as e <C 1. The solid curves denote the relation 
betwe en M and p c fo r lo 2 = 0, 0.005 and 0.01, and the dotted curve for the sequence of turning points defined by 
Eqs. d3.17| ) and ( 3.18Q . We also plot the curves of constant J where J = JK~ n (the dashed lines). According to the 
turning point theorem |^0| , a change of the sign of dM/dp c along a curve of constant value of J indicates the change 
of (secular) stability. Along curves of constant J, there are in general two turning points for different value of u> at 
p c = pi and p 2 {> Pi) for J ^ 420. The turning point at p c — p\ exists only for oj > Lu cr i t and that of p2 for u> < w cr it. 
Here, numerical calculation yields the approximate relation 



^crit - (0-35 - 



0.38)e, for 10" 4 < e < 5 x 10 



1-3 



(3.21) 



We can determine the stability for stars along a curve of constant J in the following manner. All the curves of constant 
J in Fig. 1 asymptotically approach the unstable branch of isolated spherical stars (the solid curve with uj 2 = 0). 
Stars on this branch are unstable against radial collapse. This implies that a binary star with p c > p2 is secularly 
unstable against quasi-radial collapse. Since the stability changes at each turning point, a star with p\ < p c < p2 is 
stable and one with p c < pi is secularly unstable. 
The compactness of star at X — X £ IS evaluated as 



M _ 3ek 2 
"R ~ 2(1 - 3e)/c 4 



iPc 



1/3 



F = C c eF, 



where 



kuj 2 ( n + kuj 4 / 3 \ 
F^l-—[ 1 +e) <1, 

and where C c is a constant depending on n (cf. Table I). For n — > 3, we find 

M 

— ~ 0.445eF. 
K 



(3.22) 



(3.23) 



(3.24) 



Thus, for uj — > (F — > 1; infinite separation), we recover the well-known result M/R = 0.445e for onset of the 
instability against radial collapse |13||. We note that F is less than unity for finite value of uj. Thus, the compactness 



7 



of marginally stable stars against quasi-radial collapse in a corotating binary is smaller than that for spherical stars. 
This is mainly because the stellar radius is enlarged by the centrifugal force induced by the spin angular momentum 
S. Figure 1 also shows that the instability to quasi-radial collapse sets in for a slightly smaller central density and a 
slightly larger mass than those for spherical stars in isolation, for given entropy. 

The instability for a star with p c < p\ is not due to the first post-Newtonian effect; namely, it is not associated with 
quasi-radial collapse. First, we derive the condition for the existence of such an instability. For stable equilibrium, 



the condition, d 2 E / dx 2 \M,j.K > 0, has to be satisfied. Equation ( 3.1 5| ) obviously implies that the following condition 



is necessary, independent of M and x and irrespective of the first post-Newtonian term: 



,4/3 



rj — 3kuj 4 / 3 



> 0. 



(3.25) 



Note that this condition is sufficient to ins ure st ability in the Newtonian case. 
Assuming lu <C (t]/k) 3 ^ 4 and e -C 1, Eq. (3.25) can be approximately written as 



^ < W min> 



7 1 + K(ek 2 /K) 2 / 3 Y 1 



rj — 'in(ek2/ k) 2 / 3 



1.5e(l-2.8e 2/3 ) 



(3.26) 
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V 3k 2 ep c 



1/3 



?? + /<efc 2 /«) 2 / 3 ] 1/3 

7] — 3k(£/C2/k) 2 / 3 



i?~4.6 



0.001 



1/3 



(1 + 0.9e 2/3 )i?, 



(3.27) 



where to e valua te k 2 , k, and p c /p m , we use the value for an n = 3 Newtonian spherical polytrope (cf. Table I). 
Equation ( 3.26 ) im plies that for a very small value of e (i.e., e < 0.01) with which w 2 lin < 0.01, close orbits with 



, 2 > 



see Eq. ( |3.16| )] are secularly unstable; namely, the existence of close corotating binaries may be prohibited 
when e is too small. It is noteworthy that all corotating binary stars with e — > are secularly unstable. 

In Figs. 2 (a) and (c), we show the relations for J and Q(= K n / 2 Q) as functions of lu 2 for e = 0.003 along a fixed 
value of M (here we choose M = 4.05 and 3.90). We also show the relation between J and near the minimum of J 
in Figs. 2 (b) and (d). As expected from Fig. 1, the curve of J has a minimum at lu = w m j n . For higher mas s cas e 
M = 4.05, the GR effect is important so that w 2 lin is significantly less than the Newtonian value shown in Eq. ( 3.26 ), 
i.e., ~ 0.0042. However, for M = 3.90, it approaches this value. 

Recall that the turning point theorem specifies nothing about the type of unstable mode or the growth timescale 
for the secular instability. Thus, we can suggest plausible outcomes. One of the most likely outcomes is that, given 
some dissipation, a binary with to > cj m i n could become secularly unstable to merger as in the case of stiff equation of 
state with r > 2 |fT8| . [Namely, a corotating supermassive binary appears to have an innermost stable circular orbit 
(ISCO).] It should be noted that in a binary of supermassive stars, the system evolves as a result of thermal radiation 
on the time scale r CV oi- This implies that the secular instability can play a role only if the secular timescale is as small 
as T cvo i. As discussed in Sec. II, the dissipation timescale due to gravitational waves tgw is much longer than r cvo i. 
Hence, this secular instability is irrelevant for evolution of supermassive binary stars. However, viscous damping t v ; s 
may be important. 

The minimum of J (i.e., ISCO) along a sequence of constant M appears before two stars come into contact only for 
restricted values of n (i.e., n < 1 or \n — 3| <C 1). The reason can be explained in the following manner. The angular 
momentum of one star [cf. Eq.(3.12)] is written in the form 



J(r) 



M 3/2 r l/2 



2V2 



2V2K n M 3 / 2 R 2 

5r 3/2 ' 



(3.28) 



where K n = I/(0AMR 2 ) = (5K/2)(4 7 rp m /3 / 9 c ) 2 / 3 
polytropes, n = 1, K\ ~ 0.653, and K3 ~ 0.188 fll 
of r, which is calculated as 



which is a monotonically decreasing function of n ( for Newtonian 
The ISCO is located by finding zero point of d J/ dr as a function 



dJ _ M 3 / 2 
dr ~ 2 5 / 2 rV2 



6kt. 
5 



2R 

r 



1 - 



4r dR 
3R~d7 



(3.29) 



Here, we note that dR/dr is negative, since the stellar radius is enlarg ed by the centrifugal force associated with the 
spin angular momentum which increases with decreasing r. Equation (3.29) implies that the zero point of dJ/dr for 
corotating binaries is produced by the effect of spin angular momentum and exists for r ^ 2R if K n is large enough 
> 0.8 or if \dR/dr\ is large enough, i.e., of O(Rfr), in the case that n n is small. For soft equations of state with n ~ 3, 
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K n is small ~ 0.188, but R rapidly increases with decreasing r for close orbits. (Recall that for n ~ 3, the equatorial 
radius increases significantly even with a small magnitude of the spin angular momentum; e.g., for the Kepler spin 
velocity, the equator increases by a factor of 1.5 for n — 3, according to the Roche model [pi.) This rapid increase of 
R with decreasing r produces the ISCO. 



As shown in 18 2]], ISCOs also exist for sufficiently stiff equations of state with n < 1. In this case, R does not 
change significantly with decreasing r fig] . However, k„ is of O(l) for n < 1, resulting in existence of the ISCO. For 
intermediate range of n with 1 < n < 3, neither the magnitude of n n nor the rate of increase of R with decreasing r 
is sufficient to produce an ISCO. 

From Fig. 2, we find that il has a maximum at a point where lu > w m i n . The appearance of this maximum can 
be explained in the similar manner. Along a corotating binary sequence of fixed values of M and K, p c decreases 

1/2 

with increasing uj. However, the orbital frequency fl — ojp c does not have to increase. For sufficiently small w, 
the decreasing rate of p c , dp c /cLu, due to increasing spin frequency of each star is so small that Q increases with u>. 
For close orbits, on the other hand, the effect of the spin becomes significant enough to enlarge the equator and to 
decrease p c by a large factor. Then, f2 can decrease with increasing uj. 

Assuming that the timescale for secular instability, T v j s , is as short as r evo i, we can constrain the allowed region for 
existence of supermassive s tars in corotating binar y systems. For a supermassive star of mass M, e can be written as 
1.42(M/M )" 1 / 2 [cf. Eq. Q]. Thus, Eq. ( f^g ) can be rewritten to 



M 



r>AAR [WM^) ■ (3 - 30) 



On the other hand, Eq. ( 2.22 ) has to be satisfied for achieving corotation. These two conditions constrain the allowed 



regions of M and r/R according to 



. yr \ 1/6 / fir \ -H/24 



1O 6 M J ~ R~ V 10 6 M Q , 

M < 6 x 10 6 a 4/15 M o . (3.32) 

Thus, the existence of supermassive stars in corotating binary systems may be strongly restricted provided the growth 
timescale of the secular instability is as short as t cvo i. 



IV. NUMERICAL MODEL FOR COROTATING BINARY 



To demonstrate that the analytical model presented in Sec. Ill is qualitatively correct, we perform numerical 
computation for equilibrium of corotating, supermassive binary stars. To construct a numerical model, we adopt the 
conformally flat approximation of general relativity. In this approximation, we solve part of the Einstein field equations 
assuming that the spatial line element is conformally flat (see, e.g., pl] , p2| for equations). From the post-Newtonian 
point of view, this approximation is exact only up to the first post-Newtonian order. However, supermassive stars 
are not very compact, even at the onset of radial collapse, so that the higher post-Newtonian terms do not play any 
important role. Thus, the solution derived in this formulation provides an excellent approximation. 

To demonstrate that the results in Sec. Ill are qualitatively correct, we present here only one example setting 
e = 0.003. The detailed results and detailed analysis for other parameters of the numerical computation will be 
presented elsewhere. 

The computation is performed on a 121 x 121 x 121 Cartesian grid. We here use two numerical codes which were 
previously implemented to prepare initial conditions for mergers of neutron star binaries p2^a |. Computations were 
mainly performed in the code by Shibata |p2[ , while the other code |23| was used to confirm the numerical results. 
The equations are solved in one octant (x,y,z > 0). We always assign 40 grid points for the major diameter of a 
supermassive star. Comparing the maximum mass for sufficiently distant orbits with the maximum mass for isolated 
spherical stars, which can be computed with high accuracy in a one dimensional code, it is found that with this 
resolution, the mass of the system is systematically underestimated by ~ 0.8%. Hence, our numerical results contain 
an error of this order. To check the convergence of numerical solutions, we performed computations on a 91 x 91 x 91 
grid, covering the major diameter of the star by 30 grid points. In this case, the mass is systematically underestimated 
by ~ 1.5%, indicating that convergence with increasing the resolution is achieved. 

In Fig. 3, we show the dependence of M* versus p c (compare with Fig. 1). Here, M* denotes the rest mass of one 
star. The solid lines denote the relation between M* and p c for fixed values of d = (r out + r- m )/(r 

out Tin) — 1.3, 1.5, 

and 3.05 where r out and r; n denote the coordinate distance from the origin to the outer and inner edges of star along 
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the axis which connects the centers of mass of the two stars. lo 2 (= SI 2 / p c ) is 0.006, 0.0046, and 0.00068 for d = 1.3, 
1.5 and 3.05, respectively. The dashed curves show the relation of M* versus p c for constant values of J. 

As expected from the analytical results for J > 410, the curve for constant J has one minimum at p\ and one 
maximum at pi > p\\ The former is associated with the ISCO and the latter with the onset of quasi-radial collapse. 
It is found that the turning point at p\ ap proac hes to the curve of lo 2 ~ 0.0042 for the Newtonian limit (p\ — > 0), as 
indicated by the analytical study [cf. Eq. ( 3.26 )]. 

In Fig. 4, we show the relation between J and O for M* = 4.04 near the minimum of J [compare with Fig. 2 (b)]. 
It is seen that the minimum of J and the maximum of Cl exist, illustrating that the analytical results presented in 
Sec. Ill are qualitatively correct and quantitatively fairly accurate. All these results demonstrate the robustness of 
analytical modeling for corotating, supermassive binary stars. 



V. QUASI-STATIONARY EVOLUTION OF SUPERMASSIVE BINARY STARS 

In this section we determine evolutionary sequences of corotating, supermassive binary stars to the onset of instabil- 
ity or merger. As discussed in Sec. II, supermassive stars dissipate energy by thermal radiation. We assume that rest 
mass is conserved during the evolution, since spin is far below Kepler limit in corotating systems. Although a small 
amount of mass of system may be ejected due to a stellar wind in reality J5j, we neglect such effects for simplicity. 
Since the timescale of gravitational radiation reaction is too long to dissipate the angular momentum and because of 
mass conservation, we assume that the angular momentum of the system is also conserved. We also assume that the 
effective adiabatic constant T is unchanged during the evolution for simplicity. Strictly speaking, this assumption is 
not valid since (3 depends on the entropy of the radiation which decreases during the evolution. This treatment is 
adequate only when the fractional decrease in the entropy is not very large. 

In this setting, the energy E for an equilibrium star is a function of K and x, i.e., E{K,x), and the evolution 
equation for E can be written as 

E = k 1 KMx 3/n + — x = k 1 KMx z/n , (5.1) 

OX M,J,K 

where " ' " denotes the time derivative, and where {dE / 8x)\m,j.k = since we assume quasi-stationary evolution 
along an equilibrium track. From this equation, it is immediately found that K decreases with evolution since E is 
negative. 

The decrease of K implies that M and J increase. However, J/M 2 — J/M 2 remains constant. Thus, the evo- 
lutionary sequence is specified by the value of J/M 2 . In Fig. 5, we show M as a function of p c for a fixed value 
of J/M 2 (= 28,31 and 35) and e = 0.003 as an example (thick solid lines). Along these curves, M has to increase. 
This implies that if a corotating binary initially resides in the right-hand side of the dotted curve (which indicates 
the turning points), p c increases with decreasing lo. If a corotating binary resides in the left-hand side of the dotted 
curve, p c decreases with increasing lo. Since we assume that J is con serve d, il should slightly decrease (increase) with 



decreasing (increasing) lo in these evolutionary sequences [cf. Eq. (3.14)]. Therefore, if a corotating binary initially 
resides in the right-hand (left-hand) side of the dotted curve, the orbital separation r slightly increases (decreases) 
with evolution. 

As we found in Sec. Ill, supermassive stars on the left-hand side of the dotted line are secularly unstable to merger. 
Thus, if this instability grows in a timescale as short as r cvo i, the final fate is not determined by thermal radiation 
but by the secular dissipation timescale [[[ Only when the system evolves by means of thermal emission, supermassive 
stars quasi-stationarily evolve with decreasing central density and increasing lo to merger. In either case, the final 
outcome will likely be a single, rotating supermassive star with a surrounding disk as a result of merger. Such a 
rotating star will presumably be rapidly rotating, so that its subsequent evolution may be essentially the same as for 
a single isolated rotating star as discussed in [JtJ . 

In the case when the binary resides in the right-hand side of the dotted line, the evolution will proceed via one of 
several scenarios. If the viscous timescale with regard to the tidal torque is shorter than the thermal evolution timescale 
throughout the entire evolution, corotation is preserved, and the sequence terminates at the point of unstable, quasi- 
radial collapse. We refer to this scenario as scenario (A). H owev er, vis cosit y may not be strong enough to preserve 



corotating orbits for small lo, i.e., lo < Lo v i s , where [cf. Eqs. (2.21) and (3.16|)1, 



' A point of dynamical instability is likely to arise beyond the point of secular instability, but this point cannot be identified 
by the present treatment that assumes corotation (see, e.g., 0|) 
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~ 5.6 x (m\ ^'(M) 3 " > 10-V-/3 (K) " /S . (5 . 2) 

V 10 / \10 6 M Q J \10 6 M & ) v ; 



Here we use Eq. (2.8) to derive the second inequality. In this case, the orbit may depart from corotation before 
reaching the turning point. We refer to this scenario as scenario (B). 

Scenario (A) is possible only for relatively small mass since the condition w cr it > lo V [ S is necessary. This condition 
reduces to 

M >2X10-W " \ (5.3) 



R ~ \10 6 M &/ 
M < 4 x 10 6 a 4/15 M Q . (5.4) 

Thus, only for supermassive binary stars of fairly small mass and of fairly large compactness, corotation can be 
preserved until the onset of quasi-radial instability. 

Assuming that corotation is preserved through the entire evolution [i.e., scenario (A)], we estimate the non- 
dimensional, spin parameter of each star, q = S/M 2 , at the onset of quasi-radial collapse in order to assess the 
final remnant following the collapse. Since the star in corotation has a spin angular velocity equal to fl, q is calculated 
according to 

m ( 3 Pc y 1/6 fR^ 1/2 



q = lP =K {^) [m) w (5 - 5) 



Using Eq. ( 3.22|) , we find that q at the onset of quasi-radial collapse for e <C 1 is 

q = q uje- 1/2 F- 1/2 , (5.6) 

where 

_ ^ 2fc 4 (l-3e) V /2 
qo = K { 3fc 2 ) 5 (5J) 

values of qo are at most ~ 0.4 (cf. Table I). As mentioned in Sec. Ill, u> < Lu cl -i t at a turning point for quasi-radial 



collapse (e.g., the thick dotted line in Fig. 5). With this constraint and Eq. (3.21), Eq. (5.6) reduces to 

q < 0.6q a F- 1/2 ~ 0.25F~ 1/2 . (5.8) 



Since u> is less than cj cr it [cf. Eq. (3.21)] at the onset of quasi-radial collapse, F > 0.5 [cf. Eq. ( 3.23| )]. Therefore 



if quasi-radial instability sets in for supermassive stars in corotating binaries, q has to be much smaller than unity. 
Consequently, the resulting spin is below the Kerr limit (q = 1) and the final fate is likely to be a slowly rotating black 
hole. This result is due to the spin angular momentum being effectively transferred to orbital angular momentum via 
viscosity during the evolution. 

For scenario (B), in which the corotation is not preserved for w < w v ; s , we can consider two scenarios for the 
subsequent evolution. In one scenario, each star reaches a mass-shedding limit as a result of the decrease of stellar 
radius. It then evolves along a mass-shedding sequence ejecting mass from the equator until the onset of quasi-radial 
instability 0. We refer to this scenario as scenario (Ba). As discussed in [Q, the final fate of each star after the 
onset of the instability is not clear, since q is nearly equal to unity at the onset of the instability, and rotation may 
inhibit direct collapse to a black hole. In the other scenario, each star in the binary shrinks, conserving spin angular 
momentum through the entire evolution, until quasi-radial collapse. This is the case in which the stellar radius at 
lo = cj v i s is sufficiently small to reach the turning points before reaching a mass-shedding limit. We refer to this 
scenario as scenario (Bb). In this case, q parameter could be smaller than unity so that the final outcome is likely to 
be a black hole of moderately large spin. 



We derive the criterion for the stellar radius at w = uj v i s to follow scenario (Bb). From Eq. (2.21), the angular 
velocity of each star at lu — uj V i S is 

2 4 1/2 M fM/RoY 5/4 f M \ 3/4 , x 

where R is the stellar radius at u = uj v [ s . With decreasing stellar radius R < R , uj becomes smaller than w vis . For 
u> < Wvis, the spin angular momentum of each star remains constant until the star reaches a mass-shedding limit at 
which the angular velocity approximately becomes (^-0] 
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(5.10) 



Because of the conservation of the spin angular momentum and mass, the relation, R 2 fl = Rq£Iq, approximately 
holds. Hence, 

~ 7 x 10-*a-V»^ ( (") . (5.11) 
For the second scenario, f2 has to be less than r2 s hcd at the onset of quasi-radial collapse, i.e., 

_>0.07a-»/»(^J (^-) . (5.12) 

Thus, i?o has to be less than ~ 10 times of stellar radius at the onset of quasi-radial collapse for M ~ 1O 6 M . 

In Fig. 6, we present a schematic diagram for the evolution of supermassive binaries of M — 10 6 A/q for a — 1. The 
horizontal and vertical axes denote M/R and r/R, respectively. Since M is conserved and r remains approximately 
constant in negligible gravitational wave emission, evolution proceeds at fixed M/r approximately (see solid and 
dashed lines). A binary with r/R < 2 (below the lower long-dashed line) should be regarded as a single star, so 
that only the region with r/R > 2 is relevant here. Slowly rotating supermassive stars with M/R > 6.3 x 10~ 4 and 
S/M 2 <C 1 are unstable against quasi-radial collapse (the right-hand side of the vertical long-dashed line). A binary is 
corotating if it resides in the region surrounded by three long-dashed lines. A corotating binary with r/R > 4.1 (above 
the dotted line) is secularly stable, so that if a corotating binary resides above (below) the dotted line, r/R increases 
(decreases) with evolution. (We note that for non-corotating binaries in which tidal, viscous effect is negligible, such 
secular instability would not exist, implying that r/R would always increase for them.) 

From these results, we find the following conclusion for M = 1O 6 M : (1) If M/R of a corotating binary is initially 
larger than ~ 2 x 10~ 4 , the corotation is preserved until each supermassive star in the binary reaches the onset of 
quasi-radial collapse [scenario (A)]. In this case, the product is a binary black hole of small spin parameter. (2) If 
M/R of a corotating binary is initially larger than ~6x 10~ 5 , each supermassive star in the binary reaches the onset 
of quasi-radial collapse before reaching the mass-shedding limit [scenario (Bb)]. In this case, the product is a binary 
black hole of moderately large spin. (3) If M/R is initially smaller than ~ 6 x 10 -5 , stable, supermassive binary stars 
in corotation can never be achieved [scenario (Ba)]. In this case, the outcome is not clear. If the supermassive stars 
directly collapse to black holes conserving the spin and mass, the product would be a binary black hole of large spin 
■I - I- 



VI. IMPLICATION FOR GRAVITATIONAL WAVE DETECTION 



In this section, we estimate the frequency and amplitude of gravitational waves from supermassive binary stars in 
stable circular orbits and from collapse to a black hole [i.e. for scenario (A) and (Bb)]. Possibilities for generating 
burst and quasi-periodic gravitational waves after the onset of quasi-radial instability in scenario (Ba) are the same 
as those discussed in |?J. 

Frequency / and amplitude h of gravitational waves from supermassive binary stars in circular orbits are calculated 

as 

fi / 1000AA 3/ V 10 6 Af Q \ 
/~3xlO- 6 Hz— — — ° , (6.1) 



— »(^)Go^)(^> 

where D denotes the distance from the source to detectors. Since r > 1000M for supermassive binary stars with 
M > 4 X 10 5 M Q , the frequency is smaller than ~3x 10~ 6 which is too low to be detected by LISA |^,|[. 

Close binaries will be secularly unstable to merger to a single rotating star. During merger, the coalesced object is 
likely to form a stable ellipsoid which could emit quasi-periodic gravitational waves. Even in this case, the frequency 



will be similar to that given in Eq. (6.1) and is still too low, since the radius of the ellipsoid is likely to be larger than 
- lOOOAf. 

Supermassive stars in a stable binary orbit likely collapse to black holes after quasi-stationary contraction but prior 
to merger. For the case when the orbital separation and compactness of each star are initially small, the black holes 



12 



are directly formed after onset of quasi-radial instability according to scenarios (A) and (Bb). As we found in Sec. 
V, supermassive stars are rotating with 0.3 < q < 1 at the onset of quasi-radial instability. Since the collapse is non- 
spherical, gravitational waves are emitted. Gravitational waves during formation of black holes are likely dominated 
by quasi-normal modes of black holes [B4| . Perturbation studies for the frequency of quadrupole quasi-normal modes 
provide pj| 

/- (0.03- O.OS)^/- 1 - (6-15) x lO-*(-^L-\ Hz, for < q < 1. (6.3) 

Thus, the frequency for M ~ 10 6 Af© is in the range in which LISA is expected to be most sensitive. The strength 
of the gravitational wave signal likely depends strongly on q and on the dynamics of collapse, so that it is difficult to 
assess without numerical simulation. An order-estimate yields [[llj 

V2 / 10 -3 Hz \ 1/2 /3000Mpc\ / M ^ 1/2 



where e denotes the ratio of the total emitted energy of gravitational waves to the gravitational mass energy M. 
Thus, as long as e > 10~ 10 , detection of gravitational waves will be possible by LISA in which the sensitivity around 
/ ~ 10 _3 Hz is ~ 10~~ 21 . A numerical relativistic work pj] indicates that e depends strongly on q and initial conditions. 



Since the simulations in 24 are carried out only for collapse of neutron stars, it is necessary to perform a simulation 
preparing a realistic initial condition and using realistic soft equations of state for supermassive stars to clarify the 
magnitude of e. However, Ref. |24| reports that as long as q > 0.2, e could be larger than 10 -6 . If this relation would 
hold for collapse of supermassive stars, future detection of gravitational waves for scenarios (Ba) and (Bb) would be 
possible. 

If supermassive stars eventually collapse to a black hole, a black hole binary will be formed and s ubseq uently 
decreases orbital separation radiating gravita tion al wav es d uring inspiraling in a timescale tqw [see Eq. ( 2.14| )1 . The 



frequency and amplitude are given by Eqs. (6.1) and (6.2), which show that supermassive binary black holes with 
< 200M(M/10 6 Mq) -3 / 2 emit gravitational waves of f~> 10 _4 Hz and can be strong sources of LISA. 



r 



VII. SUMMARY 



We have investigated equilibrium and stability of supermassive stars in two-body binary systems. We found that 
supermassive binary stars in close, corotating orbits with r/R < 4(M/10 6 M Q ) 1 / 6 are secularly unstable. We also have 
shown that the corotation cannot be achieved for distant orbits with r/R > 12(A//10 6 M Q )~ n / 24 , since the timescale 
for viscous transfer of the angular momentum is longer than the evolution timescale due to radiation of thermal 
energy. These facts suggest that existence of stable supermassive stars in corotating binaries is highly restricted. In 
particular for massive systems M > 6 x 10 6 M Q , corotation cannot be achieved. Evolution in such a massive binary 
does not experience transfer between orbital and spin angular momentum. Namely, the evolution of the velocity field 
of one star is not affected by its companion star. In this case, each supermassive star will be rapidly and rigidly 
rotating, as in the case of isolated supermassive stars. 

As discussed in H, rotating supermassive stars in isolation will evolve along a mass-shedding sequence to quasi- 
radial collapse. At the onset of the instability, the nondimcnsional spin parameter q is close to unity, suggesting 
that the collapsing star may not form a black hole in a straightforward manner. During the collapse, the centrifugal 
force may inhibit the collapse to a black hole near the equatorial plane to form a disk. After the formation of disk, 
non-axisymmetric instabilities may play an important role in forming a bar, or in fissioning to several blobs. On the 
other hand, if supermassive stars collapse directly to black holes, the resulting black holes will have a large spin near 
Kerr limit q ~ 1. Only numerical simulation in full general relativity can distinguish these fates. 

For supermassive stars in binaries of fairly small mass < 4 x 10 6 M Q , small initial orbital separation r ~ 4R, 
and initially small compactness, corotating orbits persist until the onset of quasi-radial collapse. In this case, the q 
parameter of each star at the onset of the instability is much less than unity since spin angular momentum can be 
transferred to orbital angular momentum by viscosity during the evolution. The final fate of the system is likely to 
be a binary black hole of small spin, although it is achieved for restricted initial conditions. 
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APPENDIX A: ROCHE BINARY 



In this appendix, we briefly investigate equilibrium and stability of supermassive stars in the Roche binary. In this 
problem, we treat binaries composed of a point body of mass m and a s uper massive star of mass M. We again assume 
that the equation of state for the supermassive star is given by Eq. (£5J]). In this setting, the total energy of the 
system E is written in the form 



E = k 1 KMx 3,n ~ k 2 M 5/3 x - k A M 7/3 x 2 
where here is determined by 

= 



Mm 



J 2 Mm 
2 + 2(M + m) ' 



M + m 



and other notations are the same as those in Sec. III. The total angular momentum is written 



Mm 
M + m 



'0-1/3 



KX - 2 n)M 5/3 , 



(Al) 



(A2) 



(A3) 



where 7/ = m/(M + m) 1 / 3 M 2 / 3 . Taking the first derivative of E with respect to x fixing M, J and K yields a 
condition for equilibrium 



= 



dE 
dx 



M,J,K 



^h, tcQ 2 M 5 / 3 

—^-KMx 3<i - fc 2 M 5 / 3 - 2k 4 M 7 / 3 x + : , 



where we use 



an 

dx 



_ 2k. 

M,J X 3 



.ILn^/ 3 + — 

3 x 2 



(A4) 



(A5) 



Equation (A4) is completely the same as Eq. ( 3.13j ), implying that we can construct equilibria in the same manner 
as that in Sec. III. 

The second derivative of E with respect to x becomes 



d 2 E 



dx 2 



9efci 



KMx 



3e- 



M,J,K 



, 7/ o 3kM 5 / 3 Q 2 f n'x 2 + nil 4 / 3 \ 

1 - 2k 4 M 7/3 

x 4 \t]'x 2 - 3«fi 4 / 3 y 



= 3ek 2 M 5/3 x- 1 - 2(1 - 3e)k 4 M 7/3 - 



3«;0 2 M 5 / 3 / T] 



'X 2 + K&/ 3 



\rj'x 2 — 3«ri 4 / 3 



(A6) 



Thus, the stability can be investigated simply by exchanging 7/ in Eq. (3.15) to 7/, so that the qualitative feature on 
the stability is essentially the same as that shown in Sec. III. A sufficient condition for existence of a stable equilibrium 
in close orbits (d 2 E/dx 2 \M,j,K > 0) is given by 



e&2 — klo 



2 ( v' + KUJ 



4/3 



where 



uj 2 = 



T] 1 - 3ku 4 / 3 

n 2 4tt p. 



3 p c 



1 



e >0, 



Equation (A7) is approximately written as 



2 . efc 2 
ui < . 



(A7) 



(A8) 



(A9) 



Thus, we find that close orbits with r ~ R are secularly unstable for a small value of e, i.e., e < 0.05(1 + m/M), in 
Roche binary. 
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APPENDIX B: EFFECT OF THE TIDAL ENERGY 



In this section, we briefly review the effect of tidal energy to the equilibrium and stability of supermassive stars in 
binary systems. In the following, we consider supermassive binary stars of equal mass and of no spin for simplicity. 
Including the tidal energy but neglecting the effect of the spin, the energy of one star is written as |l9| 

E = hKMx 3 / n - fc 2 M 5 / 3 x - k 4 M 7 / 3 x 2 + —r 2 n 2 - — - ^-=-, (Bl) 

8 2r r°X° 

where f2 is determined from dE/dr\u,j,K = as 



2M 24AM 8 / 3 . 
"=VT5-+ -a_ B ■ ( B 2) 



75 2 A ™\ /" 3/3, ■ 1 " 



and 

'■ 16 K 5j^47rp m J " (B3) 

Our notation is the same as that in Sec. Ill, except that the angular momentum of one star is written as 

Mr 2 

j = — n. (B4) 

Taking the first derivative of E with respect to x fixing M, J and K yields a condition for equilibrium 

-2 / nr ,\ 6 



= f 

OX 



_ 3k 

M,J,K n 



iKMx^ - k 2 M^ - 2 h M^x + ^! f " f , (B5 ) 

64 \&irp m ) \ r J 



where we use 



Or 
dx 



M.J 



60AM 5 / 3 A 48AM 5 / 3 \ 

1 ZjTjT- ■ ( B6 ) 



r 4™6 



For n — 3 , A ~ 0.76, and the fourth term of Eq. (B5) is written for n = 3 as 3.5 x lCP 4 M 5//3 (2i?/r) 6 . Comparing with 
Eq. (3.13), it is found that the magnitude of the tidal energy is by a factor ~ 0.02(2i?/r) 3 smaller than the magnitude 
of the spin effect for corotating binaries, implying that the effect of the tidal energy for equilibria of supermassive 
stars is much less important. 

The second derivative of E with respect to x becomes 



(B7) 



51 = ^kmx^ - 2Mfv» - 3oA f; i/3 t; +i2AM ^ 

dx 2 m,j,k n r 6 x~ x 5 r 5 - 48AM 5 / 3 

i , 7/ o 15AM 5 / 3 f2R\ e / 3p c r 2 /sV + 12AM 5 / 3 e 

= 3efc 2 M 5/3 a; - 2(1 - 3e)fe 4 M 7/3 — — — -^r- 1 ft; + ~ 

V ' 32x \ r J \4xp m ) \x 5 r 5 - 48XM 5 / 3 2 

Thus, the condition for existence of stable supermassive stars is approximately written as 

(?)■ 

In determining the stability of the non-spinning binary system, the tidal energy can important for very small e as 
e < 10~ 3 and for very close orbits with r ~ 2i?. However, as long as e > 7 x 10~ 4 , it is not important at all. 
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TABLE I. Nondimensional structure constant for Newtonian spherical polytrope with n ~ 3. e denotes 1/n — 1/3. 
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FIG. 1. Analytic mass-density relation for e = 0.003. The dependence of M on p( a ~ n >/ 2n for fixed values of lu 2 is indicated 
by the solid lines. The dependence of M on p^^ n ^ 2n fo r fixed values of J is shown by the dashed lines. The dotted line 
indicates the sequence of the turning points. 
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(c) (d) 

FIG. 2. (a) J (solid curve) and Q (dotted curve) as a function of uJ 2 along a fixed value of M (= 4.05) for e = 0.003. Two 
solid circles indicate the minimum of J and maximum of H. (b) Relation between J and Q. near the minimum of J for M — 4.05 
and e = 0.003. (c) The same as (a) but for M = 3.90. (d) The same as (b) but for M = 3.90. 




FIG. 3. Numerical mass-density relation for e = 0.003. The dependence of M on p( 3 ~ n )/ 2n f or fixed values of u> 2 is indicated 
by the solid lines. The dependence of M on p^~ n ^ 2n for fixed values of J is shown by the dashed lines. The solid dots indicate 
the sequence of the turning points. 
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FIG. 4. Numerical relation between J (solid curve) and fi for a given value of M» (= 4.04) for e = 0.003. 
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FIG. 5. Evolutionary tracks along curves of fixed J/M 2 (= 28, 31 and 35, thick solid lines from top to bottom) for e = 0.003. 
The evolution must proceed toward increasing M as indicated by arrows. The thin solid lines denote fixed values of u> and the 
dashed lines show curves of fixed value of J(= 500 (upper) and 600 (lower)), used to determine the location of turning points 
(dotted line). The thin dotted line marks the threshold between increasing and decreasing orbital separation. The thick dotted 
line marks the onset of quasi-radial collapse. 
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FIG. 6. Schematic diagram for the quasi-stationary evolution of supermassive binary stars of equal mass. The solid arrows 
indicate the direction of evolution as a result of cooling via thermal radiation. The solid and dashed lines indicate possible 
evolutionary paths. 
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